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ABSTRACT 


This study takes a wide ranging look at issues involved 
in the dynamic modeling of complex, multibodied 
orbiting space systems. Capabilities and limitations of 
two major codes (DISCOS, TREETOPS) are assessed 
and possible extensions to the CONTOPS software are 
outlined. In addition recommendations are made 
concerning the direction future development should take 
in order to achieve higher fidelity, more computationally 
efficient multibody software solutions. 
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1.0 


INTRODUCTION 


Planned space systems (large space antennas, Space-Based-Radar, 
directed energy weapons) may be called upon to provide a precision in 
maneuverability and in pointing which is orders of magnitude beyond that 
possible to date. 1 While the performance requirements are becoming more 
stringent, the spacecraft themselves are becoming more complex, larger and 
more flexible and consequently are more difficult to model. 

As well, spacecraft configurations are characteristically quite diverse from 
one another. It is for this reason that "generic" multibody codes have been 
developed 2,3 to model dynamic behavior and for use in control design. 
Reference 4 provides a good general overview, from a controls viewpoint, of 
the difficulties involved in the analysis of large flexible space structures 
(LFSS). The emphasis here, however, is on fidelity of the multibody 
modeling process and on computational efficiency. The spirit of the 
investigation is to initiate a fresh look at the major codes and to not be 
dedicated to any one methodology a priori. 

The open literature is reviewed to establish the state-of-the-art in formulation 
of multibody dynamic equations of motion and in multibody solution 
methodologies. Fundamental issues relate to the system topology are 
examined. Pros and cons of two existing codes (DISCOS, TREETOPS) are 
compared and the suitability of present structural modeling are looked at. 
Effects of other influences such as fluid motions, environmental loads and 
implications of maneuvering are also discussed briefly. As a consequence of 
the literature review, several useful extensions to CONTOPS are evident 
("free-free" flex modes, deployment, foreshortening) and an outline is 
provided for their implementation. To complete the study, some general 
conclusions are reached and suggestions are made on what might be done 
to generate more efficient, higher fidelity models in the future. 
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2.0 


MULTIBODY DYNAMIC SIMULATION 


2.1 Role of Simulation, Analysis Tools 

Closed form analytic solutions for dynamic response exist only for the 
simplest of systems: mass particles, translating finite masses, rotations of 
symmetric (or at least partly symmetric) rigid bodies, etc. As system 
complexity grows due to number of bodies, types of bodies (fluid, rigid, 
flex,...), geometric irregularity, range of motions (linear, nonlinear), the only 
recourse one has is a numerical solution. The specific list of assumptions 
grows longer as one strives for an evermore generic representation and can 
seem endless with the growth in demand for higher accuracies and the 
ensuing need to more fully capture all interaction effects. 

Ultimately situations can arise for which it is not possible to generate 
simulations capable of providing the desired accuracies. Nevertheless the 
analytic exercise remains an important one. Best-available-engineering 
models can be constructed and used to gain an understanding of important 
aspects of the problem. The modeling process itself develops insight into 
much of the basic physics involved and thus, in addition, is valuable training 
for the sorting out and interpreting of any measured data. From the viewpoint 
of controls design, this expertise forms an important building block in the 
pursuit of a workable control methodology. 

2.2 Nature of the Multibody Problem 

A multibody system is made up of any collection of interconnected bodies. In 
general, these bodies will have finite dimensions and will have mass 
distributed throughout the spatial domain. The motion dynamics tends to 
become quite complex because of kinematic/load coupling introduced via the 
interconnections. Major analytic efforts to date have been concentrated on 
determining a systematic manner for computing and keeping track of the 
forces of inertia associated with the accelerating masses. Comparatively, 
much less work has gone into establishing consistent phenomenological 
force models (elastic strain energy, dissipative mechanisms, etc.) and 
representations of environmental loads as applied to multibody 
configurations. 

2.3 Derivation of Equations 

All derivation procedures rely on some combination of D'Alembert, Newton, 
Euler, Hamilton or Lagrange principles for establishing conditions for 
dynamic equilibrium. It is taken as prerequisite that the individual procedures 
result in equations which are equivalent if not identical in mathematical form. 
Consequently this topic is not dealt with further in this study. Discussion and 
evaluations of derivation techniques can be found, for example, in 
References 5-1 4. 
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2.4 Status of Multibody Code Development 

Earlier efforts in development of multibody formalisms and computer codes 
have already been reviewed or commented on in References 5-7, 15, 16. 
Reference 17 reviews developments from a user perspective. Overall 
assessments by Jerkovsky (1978) "*5 and by Covington, et. al. (1984)‘ | 6 are 
summarized in the Appendix (Tables 1 and 2, respectively). What emerges is 
evidence of an evolutionary path being followed in the search of ever more 
realistic and hence more useful formulations. The configurations modeled 
range from the cluster and chain categories shown in Figures 1(a), 1(b) 
through to the open and closed tree topologies of Figures 1 (c), 1 (d). Some 

important factors addressed include: Degree of rigidity (100% rigid? 1 ,2 

DOF vibrations?), type of state variable (displacement? velocity? momenta?), 
choice of reference coordinates (relative? inertial?), choice of reference 
points (center of mass? fixed hinge?) and character of the interconnection 
(point? finite domain? #DOF? constraints?). 



Figure 1 Classes of multibody configurations: (a) Cluster 

(b) Chain 

(c) Open-loop tree 

(d) Closed-loop 
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The question of which formalism to adopt from the basic laws of mechanics is 
also an important one, but the exact method does not appear to be a critical 
consideration since all methodologies end up providing equivalent sets of 
equations. 5.12 Suffice it to say that for multibody applications it simplifies 
the derivation process considerably to not include nonworking forces of 
constraint. 

Another key issue in the realm of multibody modeling is centered around 
methods for dealing with the "system" aspects of the topology. For example, 
is it always advantageous to work through the barycenter of augmented 
bodiesl 8 or is the direct path approach as advanced by Ho 7 more beneficial 
all round? Hughes 6 states that, while the use of augmented bodies does 
help uncouple the system translational motion from rotations, the price paid is 
a more complicated form for the rotational dynamics. Furthermore, Reference 
6 points out that although the translations/rotations are coupled in the direct 
path procedure, the constraint forces are eliminated at the outset by choosing 
interbody hinge points as reference points and more straightforward 
interpretations are possible for quantities such as the inertia coefficients. 

All in all, a consensus (References 6,8) appears to be building in favor of the 
direct path method. This, by the way, is consistent with the methodology 
employed in the TREETOPS algorithm. 6 Still other attempts at organizing 
the system kinematic, dynamic topology is the "vector- networking" of 
Reference 1 9 and the "bond-graph" strategy followed in Reference 20. Such 
procedures are well suited to efficient implementation on a digital computer. 

It has been commonly assumed that to solve the governing multibody 
equations one must ultimately face the task of numerically inverting a mass 
matrix. There is some recent development which introduces recursive 
algebra early in the formulation to avoid this situation (References 21-23). 

2.5 Evaluation of Representative Multibody Software 

Rather than duplicate any assessments carried out to date, the focus here will 
be on comparing and contrasting essential characteristics of two of the more 
contemporary multibody codes resident in the public domain. Both 
DISCOS 6 and TREETOPS 6 capture the fully coupled, large angle, 6DOF 
translational rotational dynamics of a contiguous collection of rigid, flexible 
bodies. Vibrational displacements are handled via an Assumed-Mode type 
of expansion. In DISCOS, bodies are modeled independently and 
constrained through the kinematics and Lagrange multipliers, which also 
allows for loop closure. The net result is a fairly high order system which can 
be computationally sensitive if there are a large number of constraint forces. 
Reference 24 points out a specific problem of the closed loop formulation of 
DISCOS is the requirement that at least one member in the loop be modeled 
as elastic. TREETOPS makes use of a direct path methodology, thus 
minimizing presence of constraint forces and, in addition, the order of the 
system is adjusted according to the number of degrees of freedom specified 
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by the user. To encompass closed loop topologies, the TREETOPS has 
been modified to include specific classes of equality (e.g. prescribed 
velocities) and inequality (e.g. hardstops) constraints. This version is referred 
to as CONTOPS 25 and makes use of Singular Value Decompostion 26 to 
further minimize system degrees of freedom in the presence of constraints. 

There are other differences between these two codes. DISCOS allows for 
local body spin or locally stored momentum which TREETOPS does not. In 
DISCOS, one of the bodies can be made to follow a prescribed orbital 
trajectory whereas it is not clear that TREETOPS has this capacity. Neither 
does TREETOPS embody the aerodynamic, gravitational or thermal 
environmental loading, although this development is understood to be 
underway. Also, in DISCOS, the reference points are taken to be at the 
center of mass or at "hard" points, whereas in TREETOPS only the inner 
hinge point need be fixed. Common to both algorithms is the lack of 
provision for component mass flow (which would be needed in study of tether 
problems) or provision for coupling of axial load effects with transverse 
vibrations (see Kane, et. al. 27 or Lips 23 ) and, in both cases, body 
interconnections are assumed to be point-like. More recently some effort has 
been made to modify DISCOS so that it might include foreshortening. 23 - 30 

A significant extension to this class of multibody model is work of Reference 
31 which allows for relative translation along not a rigid, straight path, but 
rather along a flexible, curved path. 

This discussion shows how far multibody simulations have progressed. On 
the European scene, the MEDYNA 32 software is also of this calibre. As well, 
a fairly sophisticated but computationally intensive code LATDYN has been 
developed by NASA Langley 33 to deal specifically with the case of 
deployment, unfolding of the joint-dominated structures. 

2.6 Symbolic Manipulation, Paraltel Processing 

Generating equations governing multibody dynamics and the development of 
the associated computer software tends to be labor intensive and error 
prone. This situation has spawned a new approach in which the computer is 
harnessed to do both tasks using symbolic manipulation. In the United 
States, for example, SD/EXACT 34 is a symbolic code dedicated to multibody 
dynamics whereas in Europe NEWEUL 33 and MESA VERDE 36 have been 
applied to multibody spacecraft dynamics. 

Another attempt aimed at improving computational efficiencies of a multibody 
simulation is the processing of dynamic events in parallel. 37 
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3.0 


ADEQUACY OF STRUCTURAL MODELS 


The question of exactly how to best represent the structural flexibility of an 
individual body is indeed a very large subject in its own right. It is 
undergoing constant and considerable revision and we will not attempt here 
to deal with all aspects in detail. Rather the intent is to leave the reader with 
a sense of what some of the main issues are at present and to indicate 
current approaches to solutions. 

3.1 Basic Mathematical Descriptions for Flexibility 

The most common manner of representing flexibility is to express the 
displacement field as superposition of an infinite number of assumed-modes. 
These modes are normally derivable from an eigenvalue description of the 
mass, elastic properties and amounts to a standing wave model. The 
approach ultimately leads to discrete equations which can be truncated, thus 
reducing system order and rendering a transient dynamic analysis tractable. 
Truncation criteria used to be a matter of selecting the first few lower 
frequency, higher energy modes (this also eliminates high frequency modes 
which may or may not be desirable). Currently, however, the literature 
contains more sophisticated approaches which take into account, for 
example, the contributions of a given mode to rigid/flex coupling (References 
38-41). These modal cost assignment techniques have been used to arrive 
at reduced models for entire systems (see References 42-44). This approach 
can be particularly significant to a control designer faced with a high order 
finite element system model. 

The question of exactly which modal set to adopt is not clear in a multibodied 
setting because of the dynamic coupling and because of the complicated and 
not always well defined boundary condition. Ideally, the modes chosen 
would be the eigenfunctions, but they are usually not easily arrived at so that 
simplified sets of dynamic situations (e.g. constant spin) are imposed as 
reviewed in Section 1.2 of Reference 45. Perhaps it is not essential to work 
with the closest approximation to the eigenfunctions. Meirovitch, et. al. 46 
conclude that for any linear gyroscopic system it is sufficient to use any set of 
. admissible functions provided the set is complete. For a function to be 
considered admissible, it must satisfy geometric boundary conditions and be 
differentiable to order P for a system of order 2P. 

Finite Element Models (FEM) (e.g. NASTRAN) are commonly used to 
generate modal data for cases having nonuniform geometry and nonuniform 
mass, stiffness distributions. The way in which this data is used in a 
multibody code such as DISCOS has been referred to as inconsistent 
according to References 47,48. Other subtleties embodied in FEM are 
discussed in Reference 49. When dealing with very large structures, the 
order of a FEM approach can become quite large. References 50-52 discuss 
the feasibility of adopting equivalent lower order configuration models for 
those situations when a structure is built up from a series of repetitive 
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subelements. Truncated discrete standing wave solution is not the only 
modeling option. No truncation is necessary if one works in the frequency 
domain as suggested by Poelart 45, 54, but this is not possible for a 
nonlinear system. A travelling wave approach is discussed in References 55, 
56. Alternately, given the uncertainty surrounding structural characteristics, 
are stochastic models (as described in Reference 57) really more 
appropriate? 

3.2 Specialized Considerations 

3.2.1 Energy Dissipation 

Historically there has been, and continues to be, considerable seat-of-the- 
pants judgment exercised when arriving at a value to be assigned the 
damping factor of a given mode. This parameter is vital to response stability 
and settling times (particularly at high frequency) and cannot be left to 
chance if high accuracy control is mandatory. Fortunately research of a 
fundamental nature is underway as outlined in Reference 58. The objective 
is to measure damping of a small material piece of a structure, relate it, via a 
viscoelastic modeling process and time domain realization, to the full sized 
structure. Evidence of the continuing effort to rethink and experimentally 
quantify damping parameters is contained in the work of References 51 ,59- 
62, for example. Some success has been had in measuring damping in-orbit 
as cited in References 63,64 for systems having solar arrays 
(Communications Technology Satellite and Orbiter-Based-SAFE experiment) 
but for some modes the damping is much different than the predicted levels. 

3.2.2 Mass Flow Deployment 

By now we have gained some appreciation as to the complexities associated 
with configuration geometries and vehicle elasticity. Analysis becomes even 
more involved during extension or retraction of rigid and/or flexible 
components. This type of deployment introduces a variable mass distribution 
(and hence variable moments of inertia) together with additional relative 
velocities and accelerations, all factors that affect transient dynamic behavior. 
Perhaps because of its inherent complexity the problem has received 
relatively little attention. In general, available investigations tend to be more 
limited in scope than those dealing with nondeploying flexible structures as 
reviewed in Section 1 .2.5 of Reference 21 . One of the more significant 
studies (Reference 65) considers deployment of a single beam. Reference 
66 investigates the characteristics of a spinning, deploying beam in-orbit 
under the influence of the gravity gradient. One of the more general 
multibody software developments to consider individual body mass flow is 
the cluster configuration investigated by Lips (Reference 21) where it is 
shown that deployment-related Coriolis loads can be significant. More 
recently, Keat, et. al (References 16,67) are developing a multibody code 
applicable to spacecraft with unfolding and with mass flow and, as 
mentioned, Housner, et. al. (References 33,68) have developed such a code 
for joint-dominant structures. 
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Later in the report a look is taken at what the implications are of extending a 
code such as TREETOPS/CONTOPS to include deployment. 

3.2.3 Reduced Dimension Models and the Foreshortening Effect 

Reference 27 asserts that a fundamental flaw exists in multibody codes such 
as DISCOS and TREETOPS (which is not the case for the clustered 
appendage configuration of Reference 21 or the model in Reference 66). 
Concern appears to be centered around the absence of the flex- 
foreshortening-induced dynamic coupling between the axial loading and 
transverse vibration. Concern over this issue was raised earlier by Lips 
(Reference 28). In all cases the problem is described for a one-dimensional 
beam-like structure - a simplified model sometimes used to represent space 
trusses or astromast types of booms. 

Here, the literature in this area is reviewed and an attempt is made to 
determine the relative significance of this type of coupling. An indication of 
the main steps needed to introduce such an effect into the TREETOPS 
formulation follows later in the report. 

As shown in Figure 2(a), axial foreshortening refers to the translational 
displacement which a mass element experiences along the axial direction ( x ) 
due to elastic deformation (v) normal to that direction. Differential line 
element theory can be used to relate the deformed element dimension ds to 
its projected undeformed length dx as seen from Figure 2(b). This concept 
dates back to at least 1964.69 Over the years it has been introduced directly 
as a working displacement, 69,70 indirectly during evaluation of spatial 
integrals, 71 > 7 2 explicitly as an added component of translational ,66,73,74 
and implicitly as a nonlinear strain displacement effect. 75 



Tl 


Figure 2 Concept of axial foreshortening showing: a) deformation (v) 
along a transverse direction together with the constrained 
translation (ufs) along the axial direction; and (b) length of 
deformed element (ds) and its shortened length (dx) as projected 
along the undeformed x axis. 
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Kane, et. al. (1987)27 introduce effects of foreshortening explicitly as an 
added translational velocity. However, in addition to the subtleties 
associated with foreshortening, one must deal with the terminology and 
methodology unique to the Kane approach. To the uninitiated this tends to 
obscure somewhat the handling of foreshortening and the coupling induced 
through it. The problem is not helped by the absence of any reference to the 
term foreshortening or to differential line element theory. The handling of this 
problem is a new application for the Kane methodology but, at least as far as 
beam foreshortening is concerned, it does not qualify as a "new" theory. 
Nevertheless, Reference 27 is correct in pointing out that the foreshortening - 
related effects have not been incorporated into the existing multibody codes. 

The question of how significant the coupling terms are is addressed, in part, 
in Reference 76, where an approximate result is given relating the 
foreshortening - corrected frequency of a spinning beam with that of a non- 
spinning beam. The results are plotted in Figure 3 where it is evident that 
foreshortening-induced effects can be expected to become significant if spin 
rates are comparable or much higher than that frequency characterizing the 
non-spinning structure. 



NONDIMENSIONAL SPIN RATE 

(^sM) 


Figure 3 Foreshortening-induced spin stiffening effect on a 
uniform beam 
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3.2.4 


Joint Idealization 


Large flexible frame structures in space will be characterized by large 
numbers of joints each of which can present the analyst with a host of 
unknowns. Although idealized as points, joints in fact involve finite physical 
domains (a factor also in the setup of actuators, sensors). There can be 
sloppiness in joints, different friction regimes, considerable nonlinear 
behavior and flexibility. Such features could make the cost of computing the 
response prohibitive. Only recently has attention been focused on this area. 
Perhaps the most advanced effort to date is the Large Angle Transient 
Dynamics (LATDYN) software under development at NASA Langley 
(References 33,68). Independent work of Reference 77 suggests that the 
presence of joints may not alter modal frequencies significantly provided the 
overall structure is not too short. Some work is also being done to identify 
energy dissipation, taking place at the joint (see References 78-80). 
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4.0 


OTHER FACTORS IMPACTING FIDELITY OF MULTIBODY 
SIMULATION 

Fluid Interaction 


4.1 


Just as elastic systems make it difficult to quantify response, the presence of 
fluids can present an even more vexing modeling problem being particularly 
sensitive to temperature, pressure. With the advent of liquid apogee motors, 
the liquid/solid ratio can now be as high as 55%, and hence the problem is 
commanding much more serious attention both in Europe and in the U.S.A., 
not to mention SDI applications. An overview of the state-of-the-art can be 
had by reviewing the papers given at the colloquium dedicated to this 
problem and cited in Reference 81. 

A number of presentations concentrated on development of fundamental 
models applicable to rotating, free-surface, fluid behavior. However, due to 
the complex nature of the problem only restricted approximate solutions, in 
the form of one homogeneous eigenmode for example 813 are forthcoming. 
As well as analysis, a significant amount of the work is dedicated to 
experimental testing of ground-based gyro-fluid systems (e.g. Vanyo, 
Reference 81 b). While the test data appears to be reliable, it is not certain 
that fluid character observed at ground level will be the same in orbit. To 
address this question, ESA has observed fluid behavior in orbit on the 
Shuttle SPACELAB 1 mission 81 c . 

An experimental method which appears to give good agreement (better than 
20%) with in orbit measurement, is the ground-based "Drop-Test" 81 d . A 
dynamically scaled model of the satellite/fluid system is allowed to free-fall. 
Nutation data is transmitted continuously throughout the duration of the drop. 

In addition to generating analytic fluid models and measuring fluid/system 
behavior, serious efforts have been made to construct overall system models 
81a, 81 e. McIntyre of Hughes Aircraft outlines an approach that is analogous 
to that used when dealing with flexibility effects. Attitude equations are 
formulated for the spacecraft which couple directly to fluid motions which, in 
turn, are characterized in terms of discrete coordinates (derived from a 
truncated Assumed-Mode expansion of the fluid displacement). Required is 
a set of representative fluid mode shapes which can be arrived at either 
analytically, experimentally or in combination. Such a method, when 
generalized, can also model fluid/flexibility interactions (platform asymmetry, 
etc.). Another effect which can become large as fluid mass ratios increase, 
and which can effect the system dynamics, is a shift in center of mass 
location. 

Attempts to come to grips with this complex phenomenon are continuing at 
both the experimental and theoretical levels as evident from References 82, 
83 which deal with coning cylindrical fluids and with clusters of gaseous 
bubbles contained within an oscillating column of liquid, respectively. 
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4.2 


Environmental Loads 


An analyst may go to great lengths to generate a high fidelity model of the 
inertia-related forces of a vibrating gyroscopic system only to find that 
external forces applied by the environment are known only very 
approximately in some instances. References 84, 85 are among the more 
recent efforts aimed at improving atmospheric density modeling, whereas 
Reference 86 is one of the latest attempts to come to grips with solar- 
radiation-pressure-induced loads. To develop the best possible 
environmental modeling, one should also consult with the comprehensive 
bibliography on the subject contained in Reference 87. To date, generic 
multibody environmental models are not available since the description used 
is linked closely to the kinematic gendre adopted within a given formalism. 

4.3 Requirement for Controlled Maneuvers 

It is expected that upcoming space systems will be required to exhibit a 
higher degree of maneuverability. For example, the entire spacecraft (or part 
of it) may have to be re-oriented rapidly or it may have to be spun up or spun 
down in the shortest time while keeping vibrations to a minimum. Reference 
88 discusses open loop strategies for vibration control of an Orbiter-based 
experiment having long (up to 150m), highly flexible booms (Note that in this 
case flex foreshortening is included at the outset.). References 89-92 also 
point out strategies for slewing flexible systems so as to minimize the post- 
maneuver oscillation. References 93-94, on the other hand, deal with 
optimal attitude control for large angle maneuvering of rigid vehicles. None 
of these studies yet rely on any of the major multibody codes. 

4.4 Numerical Considerations in the Solution of Mixed Differential- 
Algebraic Equations 

Several methods have been developed to solve sets of mixed Differential- 
Algebraic Equations, hereafter abbreviated DAE. Gear 95,96 has shown 
that his stiff numerical integration algorithm can be used to solve certain 
kinds of DAE. Certain other DAE can be converted into equivalent forms that 
are solvable by Gear’s method. There exist, however, systems of DAE that 
cannot be solved by this method. In References 97-100, Gear's idea has 
been successfully employed for analysis and optimal design of mechanical 
systems and electrical networks. Mathematical models of dynamic 
mechanical systems, using Gear's method of solving DAE, make the 
equations of motion stiff, thus involving considerable computational cost and 
time. Furthermore, this formulation requires several problem dependent 
derivative expressions that make automatic computer generation of the 
equations of motion difficult 98,99 | n fact, Reference 101 shows that DAE 
systems are not ordinary differential equations at all. 

A second approach, due to Baumgarte, 102 uses ideas from feedback control 
theory to construct a modifed differential equation that implicitly accounts for 
the constraint equations. A difficulty, however, is the selection of factors in 
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this method that greatly influence accuracy of solutions. Proper values of 
these factors are not known and experts are in disagreement over suitable 
values. Furthermore, there is no way to impose positive error control on 
constraint violations for this method. Work is in progress to develop methods 
to determine these factors, based on constraint violations that ensure that the 
constraint violations at the next time step are minimized. The major 
drawback of these methods is that all n generalized coordinates must be 
integrated, whereas only n-m of them are independent. Furthermore, m is 
often large and n-m is small. Thus, it would be attractive to integrate for only 
n-m independent generalized coordinates. Solution for the remaining m 
dependent generalized coordinates from constraint equations can then 
ensure positive error control on constraint violations. The n-m generalized 
coordinates that are computed by integrating the equations of motion are 
called independent generalized coordinates and the remaining generalized 
coordinates are called dependent generalized coordinates. 

Independent generalized coordinates can be picked from the set of n 
generalized coordinates that describe the system configuration. Wehage 
and Haug 103,104 developed an algorithm to pick an acceptable choice of 
independent coordinates, using LU decomposition of the constraint Jacobian 
matrix to identify independent and dependent coordinates. This algorithm 
performs satisfactorily, but occasionally leads to poorly conditioned matrices. 
It is then required that LU factorization be repeated and a new set of 
independent generalized coordinates be selected. The result is an increase 
in computing time and greater propagation of integration error than desired. 

In Reference 26 a method based on singular value decomposition is used to 
solve the equations of motion of dynamic systems. This method, however, 
does not provide physical insight into the solution process and is not 
convenient for automatic generation of the equations of motion for general 
mechanical systems in a general purpose computer program. 

A broader class of vectors p of independent generalized coordinates can be 

p = VI x q 
(n-m)xl (n-m)xn nxl 

where VI is a transformation matrix. In the algorithm of References 103,104 
matrix VI is a Boolean matrix. Therefore, only individual generalized 
coordinates are chosen as independent. Use of a Boolean transformation 
matrix to pick independent generalized coordinates from q is rather 
restrictive. It is desirable to enlarge the set of possible independent 
generalized coordinates by allowing matrix VI to be a more general, possibly 
nonsparse matrix. Such a choice of VI can allow linear combinations of 
physical generalized coordinates to be selected as independent. If the rows 
of VI are mutually orthogonal, the resulting linear combinations of individual 
generalized coordinates will be mutually independent. Specific forms of 
such a transformation can be chosen, based on properties of the DAE under 
consideration and desired properties of the transformed coordinates. 
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5.0 


EXTENSIONS TO THE CONTOPS ALGORITHM 


This section discusses revisions that, if fully implemented, could increase 
both range of applicability and efficiency of the CONTOPS code. 


5.1 Generalization of the CONTOPS Hinge 

Currently, there are two restrictions imposed in CONTOPS. 

(i) The reference point hk (Figure 4) of the reference body is its mass 
center, 



Figure 4 Geometry of deformable body j. 


(ii) Modal displacements of reference point hk and modal rotations of frame 
bk (Figure 4) are assumed to be zero for all bodies. For the reference body 
the second restriction is naturally satisfied if the free-free mode shapes are 
utilized in defining its elastic deformations. For all other bodies in the system 
the second restriction is naturally satisfied by the fixed-free modes of the 
individual body (fixed at point hk in Figure 4). CONTOPS results may not 
yield desired accuracy if the modeshapes are obtained with boundary 
conditions other than those described above. 

Because the mass center is taken as the reference point for the reference 
body two vector quantities Jiki and )k defined by equations 18,22 of 
Reference 3 (see attached copy) have been assumed null (for free-free 
modes of the reference body this is a valid assumption). With this exception 
the form of the equations is same for all bodies in the system. The first 
restriction may be removed by inserting non-null vectors for Jikj and )k for the 
reference body Bk . 
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Removal of the second restriction is not an easy matter. It is pointed out in 
Reference 3 that the elastic deformation, u k . is measured with respect to 
frame bk . Furthermore, it is required that the elastic deformations remain 
small. Thus the frame bk and point hk cannot be arbitrarily defined. In what 
follows the formulation of Reference 3 will be modified which will allow a 
CONTOPS user to: 

(i) define reference point and reference frame for individual bodies 
independent of the definition of hinge. 

(ii) generalize the existing hinge definition to accommodate a wider class of 
boundary conditions for the mode descriptions of the substructures. 

Figure 5 describes the kinematics of body Bk and the kth hinge. Frames pk 
and qk and vector yk fully define the kth hinge. The reference point ck for Bk 
is arbitrarily chosen. 



Figure 5 k 1 * 1 hinge and k 1 * 1 body kinematics 
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Following the notations of Reference 3 from Figure 5. 
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Equations 3,4,5 are sufficient for identifying the coefficients of generalized 
speed (equations 26,27,28 of Reference 3). 


Figure 6 here would replace Figure 4 of Reference 3. 



Figure 6 Vectors defining coefficients of generalized speed. 


Replace equation 26 of Reference 3 by 
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Replace equation 27 of Reference 3 by 
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(7) 


Replace equation 28 of Reference 3 by 
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Equations 14,15 of Reference 3 remain intact except superscript hk is 
replaced by ck . 

5.2 Implications of Deployment and Flex Foreshortening 

As currently formulated, the CONTOPS multibody algorithm does not allow 
for mass deployment or retraction of individual members nor does it take into 
account the geometric foreshortening that accompanies flex deformation. In 
light of this, the existing recursive equations are reviewed once again. This 
report discusses some basic considerations involved in modifying CONTOPS 
to incorporate such effects. 
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5.2.1 


Existing CONTOPS Kinematics (Reference 3} 


The kth body (Bk), shown earlier in Figure 4, is in the deformed state. 
Elemental mass dm in Bk is located relative to an inertial reference by the 
vector (using Equation 1): 



R K locates the hinge point relative to inertial 
r k is fixed in B k 


0 ) 



represents 


elastic deformation 


As already implied, in order to arrive at a solution in terms of a set of ordinary 
differential equations (and to pave the way for truncation to an approximate 

u k = EiJ (r*) , £t) (10) 

1 

Overall inertial velocity of dm then is given by 


R k = tfk + „ k x (r k + u k ) + u k (11) 


where 

( ) = inertial time derivative. 

o 

( ) = differentiation in Bk - fixed frame. 
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5.2.2 


Modifications due to Changing Mass 


It is common in satellite design to store elements in as compact a manner as 
possible prior to flights and then to unfurl the device once in orbit, allowing it 
to extend to full size. Such Storable Extendable Members (STEM) 
complicate dynamic studies by introducing a time variable mass feature into 
the analysis which becomes even more complex when components are 
flexible. To allow for this in the CONTOPS type of multibody code implies the 
following (Reference 21,28): 

1. A structural mass element can be viewed asbeing convected by a 
prescribed rate of deployment (Y D ) in a manner analogous to a fluid 
continuum model. Then any local time derivative must be augmented 
so that 

d( ) | = iU + ^d‘ v) ( } 

d t I kth at no spin (12) 

body no deployment 

2. When differentiating the discrete expansion given by Equation 10, it 
must be recognized that u K is not zero but now is implicitly time 
dependent as a result of the convective motions discussed above. 

3. Mass must be conserved at all times. This leads to another equation 
which governs the motion - that is, the equation of continuity. 

4. Mass flow changes system properties such as moment of inertia, 
momentum. 

Implementation of these considerations is simplified if, rau.er than a general 
continuum, one analyzes a one-dimensional beam-like substructure, for 
example. If it is uniform and if only linear vibrations are involved, then 
deployment velocity is identically equal to rate of change in length as 
measured along the undeformed axial direction. 
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5.2.3 Incorporating Flex Foreshortening 


5.3 


To include foreshortening (ufs) in a CONTOPS type of algorithm, it can be 
added explicitly to the displacement field 21 so that Equation (9) becomes 


R k = 


R hk + r k + u k + u k 

— — — IS 


(13) 


To date ufs has been approximated only for beam-like, plate-like structures. 
If we take, 




]c 

+ U 2 —2 
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(14) 


where unit vector ai lies along the undeformed axial or 'long' direction of the 
structure, then 


k' 2 


Ufs * [1/2 J* [(u* ) 

where ( )’ = 3( ) /dx * 


(u k ') 2 ] 


dx ... ] a. 


(15) 


+ higher order terms 


When differentiating to arrive at velocity, acceleration it must be kept in mind 
that the approximated foreshortening displacement is time dependent. 
Alternatively, since the TREETOPS approach makes use of the partial 
velocity description in the kinematics, one can adopt the procedure given by 
Kane et. al 27. in this case the foreshortening is embedded in the velocity 
field. 


A Model Reduction Approach for CONTOPS 

The CONTOPS formulation is based on the method of substructures. The 
elastic deformations are expanded in terms of the natural modes of individual 
flexible substructures. The number and type of substructure modes required 
to accurately describe the elastic behavior has not always been addressed 
properly. Unwise choice of mode shapes will result in modeling error 
producing poor quality solutions. Also, retaining the substructure modes with 
insignificant contribution to the solution results in unnecessary computational 
penalty. 

In what follows we outline an approach based on the modal cost analysis 
(MCA) approach, a well proven technique for model reduction. The MCA 
approach to model reduction is widely used in the control science discipline. 

This approach was first proposed by Dr. R.E. Skelton of Purdue University 
(See Ref. 38, for example). In this approach a system is defined to be 
comprised of many components. The central idea of the modal reduction 
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method of the MCA is to exploit the precise statement of the optimal control 
problem in order to predict which system components will make the largest 
contributions in the total quadratic performance criterion. These components 
are retained and the balance are discarded. Mathematically, the concepts of 
this model reduction process are explained as follows. Let a state space 
description of a time invariant system be given as 

x(t) = Ax(t) + Dw(t) 

y(t) = Cx(t) ( 16 ) 

where the n states are represented by x(t), the k outputs are y(t) and the q 
process noise are w(t). w(t) is assumed to be zero white noise with unit 
intensity. It is assumed that the model is controllable, observable and 
asymptotically stable. It is also assumed that there are k independent 
outputs, i.e., rank (C) = k and q independent inputs, i.e., rank (D) = q. 

Let a quadratic cost function associated with ( 16 ) be given as 

V = lin E[V( t ) J , V(t) - \ f t y(T) T Qy(t)dT (17) 

t-)« *• Jq 


where E is the expectation operator and Q is a symmetric positive definite 
weighting matrix. 

Now, partitioning the state vector x(t) = (xi (t), X2(t) x n (t)) T , the component 

costs Vj associated with each component Xj is defined by 


1 11 . 


1 = 1 , 2 . 


( 18 ) 


and are calculated according to the following component cost formula 
V A “ [C T QCXJ 11 . i - 1. 2 

where (.)jj denotes the (i,i) element of (.) and where X, the steady state 
covariance of the states, is defined by 

I “ li« Efx(t)x(t) T ] 

t->® 


(19) 


( 20 ) 


satisfies 
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0 


( 21 ) 


ia t + ax + dd t = 

Clearly, since V = Tr(CTQCX) the component costs satisfy the cost 
decomposition property 

n 

V = J Vj (22) 

i=J 


The component costs are ranked in the following manner 

v l ^ v 2 > . . . > v n (23) 

and the state vector is reordered correspondingly. The most critical 
component of the system is xi having component cost V-| and the least 
critical component of the system is x n having component cost V n . 

Based on this component cost ordering, components having lower 
component cost are truncated. 

The MCA approach outlined above may be applied to CONTOPS for 
selecting the appropriate modes from a given set of modes for individual 
substructures. 

For CONTOPS application the state vector x(t) may be arranged as follows 
*(t) * I.... 

ath. (a+l)th ( 24 ) 

entry entry 


wherenjttJis the 1th modal coordinate of the kth substructure (body B^). The 
component cost corresponding to the 1th mode of B[< is 


v(i^.n^) 


Ic t qcx] 


a, a 


+ [c t qcxj 


(a+1 ) , (o+l) 


(25) 


where X satisfies equation (21). In the technical discussion above the 
component cost computation is based on a scalar, y T Qy. This concept may 
be expanded to account for multiple objective problem. This can be 
accomplished by computing the component costs for each of 



i *■ i. 


k 


(26) 


The component costs for each y2-| should be ranked and the resulting 
dominant modes based on each y 2 1 should be retained. The MCA algorithm 
can easily be integrated into CONTOPS. 
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6.0 


DISCUSSION AND CONCLUSIONS 


A wide range of issues related to the modeling of the dynamics of multi- 
bodied systems has been reviewed In this report. While it has not been 
possible to evaluate all areas fully, it is hoped that bringing so many of these 
areas together in one report has served to demonstrate the scope of the 
problem. 

Casting a multibody dynamic system into an understandable, manageable 
framework is a challenging task because of the large numbers of bodies 
involved and because of the high degree of coupling between them. 
Nevertheless, the exact procedure adopted for derivation of the equation set 
does not appear to be a concern since all aproaches will, if competent, 
produce equivalent mathematical representations. For instance, equations 
can be put into a 1 st order velocity format using either a Newton-Euler 
approach or using the procedure favoured by Kane and Levinson 10. The 
degree of effort and preference for a method will, in large measure, always 
be dictated by the analysts experience. Choice of coordinates and frames of 
reference may leave more room for debate for specialized applications but 
for a very general model this is not so much the case (6DOF is 6DOF!). What 
does stand out as a predominant issue facing multibody experts is the choice 
of a configuration topology and of methods for dealing with it. In this regard it 
would appear that the direct path method (which is conducive to use of 
recursive interbody equations and makes it straight forward to eliminate 
nonworking constrained forces) is gaining acceptance over the barycentric 
approach. One of the more versatile topologies is the tree analogy with open 
or closed branch loops. The bookkeeping aspects of the kinematics might 
also be helped through vector-networking or through bond-graph models. 

Deriving equations need not be done manually. According to the literature, 
a digital computer can be used as a symbol manipulator to generate 
equations based on whatever rules of mechanics one chooses and, if 
desired, to translate these equations directly into computer code (e.g. 
FORTRAN) for solution. This approach can be expected to yield more and 
more efficient computational algorithms and to be less error prone as a 
problem grows in complexity. A decided advantage here is the capacity to 
generate equations for the specific problem at hand rather than to generate 
one algorithm which attempts to cover all possible situations. Numerical 
efficiency can be further enhanced by carrying out computations in parallel 
whenever possible rather than in a serial manner as has been more the 
custom to date. This strategy would appear to complement the order n 
recursive modeling which avoids the need to directly invert a system order 
mass matrix and hence should augment further the speed at which the 
simulation can be executed. 

Obviously the order of a model has a strong influence on computational 
performance. As alluded to already, the order n approach combined with 
symbolic manipulation can provide a much reduced model for the system. 
Two other factors are important in determining system size - flexibility and 
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constraints. Equivalent continuum modeling together with modal cost 
analysis or simple truncation can be used to minimize the flex contribution. 
Singular Value Decomposition can be applied to minimize degrees of 
freedom for constrained systems. Model reduction can then be applied to the 
resulting system. This combined reduction may yield the most meaningful 
and efficient definition of the system dynamics depending on the 
computational penalty paid to achieve each stage of the reduction. 

The algorithm adopted to carry out a numerical integration also plays an 
important role in how accurately and how quickly a particular piece of 
software can respond. While this has not been addressed in great detail 
here, it is implicit throughout. Evaluations should be made to etermine 
relative merits of Runge-Kutta versus Predictor-Corrector methods. The Gear 
approach can be vital to the success of a 'stiff' multibody situation in which 
characterise times associated with the motions are highly divergent (this can 
occur easily for cases involving extension of a substructure from a very small 
length say out to a very large value). 

It has been pointed out that it is not the fundamental principles of mechanics 
which stand in the way of achieving a high level of modeling fidelity. From 
this review it is concluded that, in a very general sense, what does restrict the 
modeling, simulation process is: 

(i) level of geometric detail used to describe a finite dimensioned physical 
system; 

(II) adequacy of the phenomenological models. {This ranges from 
uncertainties present in the basic structural model (including, in 
particular, energy dissipation) through to questions related to fluid 
behaviour and including effectiveness in accurately capturing the effect 
of environmental influences.} 

(iii) efficiency of the solution algorithm; and 

(iv) numerical techniques (integrators, parallel vs serial processing, ...) 

It is clear from this investigation that the state-of-the-art has matured to the 
point where a contemporary generic multibody code can be assumed to 
include the following features: 

• arbitrary number of point-connected bodies 

• allows large angle relative interbody rotations (3DOF) 

• allows large relative interbody translations (3DOF) 

• allows for linear vibrations (3DOF translation) 

• some form of loop closure exists 
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While DISCOS and TREETOPS/CONTOPS fall into this category, they do not 
routinely account for deployment, for example, TREETOPS is superior from a 
numerical viewpoint in that, as much as possible, constraint forces are 
eliminated. The DISCOS Lagrange multiplier feedback approach, which 
includes a large number of constraint forces, can be expected to provide for a 
much more rapid build up of the numerical error. Further, the size of system 
matrices in TREETOPS is not fixed by the number of hinge points and bodies 
(as in the case of DISCOS) but is determined only by the number of degrees 
of freedom input to the program for a specific application. DISCOS, on the 
other hand, has the advantage of permitting prescribed motion (including 
orbital) and body spin (or locally stored momentum). This demonstrates that 
the search for a completely unrestricted 'generic' code is as allusive as ever. 

One area neglected in both of the codes discussed above, and in other 
known general purpose software (except for the development in Ref. 21), is 
the effect of nonlinearity in the strain displacement relationships of a flexible 
body which gives rise to the foreshortening in the displacement field as 
discussed in Section 3.2.3. The significance of this effect is summarized 
qualitatively here by paraphrasing Ref. 28.: "The generality of a formulation 
is compromised if one is interested in simulating response for space systems 
which are large physically, which undergo large attitude excursions (in terms 
of both displacement and rate) and which have vibrating appendages that 
engage in significant local translations if the effect of interaction of axial 
loading with the transverse vibrations is not accounted for. It is well known, 
for example, that centrifugal loads associated with rotation of a beam-like 
appendage affect terms which are linear in the vibration variables." In this 
report, it has been shown quantitatively that the spin-stiffening becomes 
important as spin rates approach or exceed the fundamental frequency 
characterizing the nonspinning structure. Note that this is the same effect 
focused on most recently by Kane et al. 27 Their claim that it has not been 
included in the major multibody codes is a valid one, however, the 
phenomenon itself is not new and has already been dealt with extensively in 
the general literature. Fortunately, existing codes can be modified to 
accommodate this effect without completely rederiving equations or 
completely reprogramming the existing codes. 
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7.0 


RECOMMENDATIONS 


This investigation has given insight into the level of complexity and the level 
of uncertainty surrounding multibody simulation dynamics. A series of 
recommendations is discussed here aimed at providing for improved 
dynamic models or improved computational efficiency of either present-day 
or of next-generation software. 

The foreshortening effect forms a valid part of the kinematics and should be a 
part of existing multibody codes and should be allowed for in any future 
modeling developments. It can require considerable additional calculations 
even though it is not significant in every application. For this reason, it should 
be retained in the form of an option. Note that CONTOPS can likely be 
amended more easily through the generalized speed than can DISCOS 
which relies on kinetic energy. 

The fact that the major codes are amenable to the type of revision alluded to 
above is encouraging, since it implies that significant enhancements can be 
achieved without generating an entirely new program. Every effort, therefore, 
should be made to capitalize fully on existing codes before contemplating 
construction of any new package. On the other hand, any single piece of 
software should not be relied on to provide the best solution for all 
applications. Rather it is more realistic to expect each code to apply to a 
certain class of problems and configurations. What is important is that the 
code description spell out clearly the known capabilities as well as the 
known limitations of a given formulation. 

Regardless of which code is being worked with or whether a new model is 
being generated, in order to arrive at higher fidelity synthesized models it is 
crucial that the best possible phenomenological models (impedance, 
standing wave, traveling wave, stochastic, viscoelastic ,...) of substructures 
and components be first determined and then implemented. Of course, 
because of the potential for dynamic interactions, all known subsystems 
should be modeled as well Another very basic factor that should not be 
overlooked is the level of geometric detail used in describing the physical 
domains and interbody connections. The more complete the representation, 
the more realistic the model and the less likely it is that important coupling 
effects will be lost. This can become particularly important when dealing with 
very large frame, truss space structures. 

An area that will always be of interest is the manner in which the fundamental 
principles of classical mechanics are interpreted and applied in order to 
derive the governing equations. At this point in time this represents a fairly 
mature part of the technology. Consequently it is recommended that the 
emphasis now be shifted toward topological considerations and to questions 
relating to the efficiency of the solution algorithm itself. For example, this 
study indicates that use of a direct path method of formulation which 
eliminates constraint forces whenever possible is preferable for that based 
on use of barycenters and augmented bodies. Graph theory should be relied 
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on to provide a method of accounting for interbody connectivity and related 
kinematics. A newer approach that merits attention for any new development 
is the recursive, order n algorithm which avoids the need for a direct 
inversion of the system mass matrix. 

Another tool that should be developed further to aid the analysis is symbolic 
manipulation. This technique harnesses the power of the computer to 
provide more rapid, less error prone sets of equations for complex 
dynamically coupled systems. 

From a numerical viewpoint there are many challenging issues. In order to 
bring the system dimension down to a tractable order, model reduction 
should be applied at both the subsystem and at the system level. The order 
of constrained systems can also be minimized through techniques such as 
singular value decomposition. As with model reduction, however, there are 
computational costs which have to be factored in to see if real savings are 
achieved overall. 

The accuracy and efficiency of numerical integrators (Runge-Kutta, predictor- 
corrector,...) needs further evaluation for space dynamic applications. How 
do the methods compare in general? What are the real benefits of having a 
variable step size in the presence of controllers? A major issue is how 
effectively can the "stiff" dynamics be dealt with. Is accuracy degraded 
seriously in presence of constraint feedback? 

Computational efficiency can be enhanced by carrying out calculations in 
parallel whenever possible rather than following the step by step serial 
approach so common in the past. In addition to parallel processing, one 
should take advantage of symbolic manipulation to write the necessary 
computer code. This promises to be faster and should definitely cut down on 
debugging costs. It is likely that this can be used when extending existing 
codes. 

What is clear from this study, is that considerable uncertainty persists in the 
modeling of multibodied dynamic systems. Ideally one would like to verify 
the working of simulation software by experimental measurement. This is not 
yet practicable for the full scale in-orbit environment. Hence, it is 
recommended here that a set of space-related baseline simulation test cases 
be established for use in evaluating performance of a given code. The test 
cases should include simplified, idealized configurations (e.g. rigid dumbell) 
through to actual mission configurations (representative communications 
satellite, Orbiter and payload, spinner;...). This process does of course raise 
questions of its own. Nevertheless, it is viewed here as a very important step 
forward in providing a framework for validating complex multibody codes and 
for establishing their range of application. 
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APPENDIX 

PREVIOUS ASSESSMENTS 
OF MULTIBODY FORMULATIONS 



TABLE 1 


Jerkovsky - Reference 15 



Table 1 Comparison of 10 n-body dynamics formulations 


n-bodv configuration 

Tree of rigid bodies; 
no relative translation. 


Tree of rigid bodies; no 
relative translation (ex- 
cept for point masses). 


Cluster of rigid bodies; 
no relative translation. 


Tree of rigid bodies; no 
relative translation. 


Tree of rigid bodies 23 ; 
terminal bodies may be 
flexible; 24 hinge points 
may be time-dependent; 
no relative translation. 

Tree of rigid bodies 23 ; 
terminal bodies may be 
flexible 2 ® no relative 
translation. 

Tree of flexible bodies; 
closed loops treated via 
Lagrange multipliers 30 ; 
relative translation al- 
lowed. 

Tree of rigid bodies; 
terminal bodies may be 
flexible; no relative 
translation. 

Chain of flexible bodies; 
no relative translation. 


Arbitrary configuration 
of flexible bodies; closed 
loops allowed; relative 
translation allowed; pre- 
scribed motion allowed. 


Dynamics state variables 

Inertial angular velocities 
of each body; inertial 
linear velocity of com- 
posite. 


Relative angular veloci- 
ties between bodies; 
inertial linear velocity of 
composite. 


Relative gimbal angle 
rates between bodies; 
inertial linear velocity of 
composite. 

Free components of 
angular momentum of 
outward composites. 


Relative gimbal angle 
rates between bodies; 
inertial linear velocity 
of a material point of 
Body 1. 

Relative gimbal angle 
rates between bodies; 
inertial linear velocity of 
composite. 

Relative gimbal angle rates; 
relative displacement rates; 
inertial linear velocity of 
composite. 


Relative gimbal angle rates; 
inertial linear velocity of 
a material point of Body 1. 


Relative gimbal angle rates; 
inertial linear velocity of 
a material pint of Body 1. 


Inertial angular velocities 
of each body; inertial linear 
velocity of material point 
of each body. 


Comments 

Uses barycenters and 
augmented bodies; 
constraint torques 
obtained via Lagrange 
multipliers. 

Use of relative angular 
velocities result sin for- 
mation of composites; 
constraint torques 
removed by use of 
symmetric projectors. 

Constraint torques do 
not appear. 


Constraint torques do 
not appear, mass matrix 
is not computed 
explicitly. 

Equations are obtained 
inductively; constraint 
torques do not appear; 
uses a nonsymmetric 
"mass matrix." 

Uses barycenters and 
augmented bodies; con- 
straint torques do 
not appear. 

Uses barycenters and 
augmented bodies; 
constraint torques do 
not appear except with 
closed loops. 

Does not use bary- 
centers and augmented 
bodies; constraint 
torques do not appear. 

Uses quasistatic modes 
plus vibration modes; 
constraint forces and 
torques do not appear. 

Free body equations are 
written for each body; 
all constraint forces and 
torques are obtained via 
Lagrange multipliers. 


Date.Authors 

'65 Hooker, Margulies 
'67 Roberson, Wittenbui-g- 


'67 Velman 


'67 Palmer 


'69 Russell 


'69 Farrenkopf 
'71 Ness, Farrenkopf 


70 Hooker 
73 Likins 


72 Roberson 
74 Wittenburg 
74, 75 Boland, Samin, 
Willems 


74, 75 Frisch 
74, 77 Ho 
75 Hooker 


74 Ho, Hooker, 
Margulies, Winarske 


75 Bodley, Devers, 
Park 



